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ABSTRACT 

Context. Contrary to supermassive and stellar-mass black holes (SBHs), the existence of intermediate- mass black holes 
(IMBHs) with masses ranging between 10'^"'' Mq has not yet been confirmed. The main problem in the detection is 
that the innermost stellar kinematics of globular clusters (GCs), the natural loci to IMBHs, are very difhcult to resolve. 
However, if IMBHs reside in the centre of GCs, a possibility is that they interact dynamically with their environment. 
A binary formed with the IMBH and a compact object of the GC would naturally lead to a prominent source of 
gravitational radiation, detectable with future observatories. 

Aims. We use A'^-body simulations to study the evolution of GCs containing an IMBH and calculate the gravitational 
radiation emitted from dynamically formed IMBH-SBH binaries and the possibility that the IMBH escapes the GC 
after an IMBH-SBH merger. 

Methods. We run for the first time direct-summation integrations of GCs with an IMBH including the dynamical 
evolution of the IMBH with the stellar system and relativistic effects, such as energy loss in gravitational waves (GWs) 
and periapsis shift, and gravitational recoil. 

Results. We find in one of our models an intermediate- mass ratio inspiral (IMRI), which leads to a merger with a 
recoiling velocity higher than the escape velocity of the GC. The GWs emitted fall in the range of frequencies that a 
LISA-like observatory could detect, like the European eLISA or in mission options considered in the recent preliminary 
mission study conducted in China. The merger has an impact on the global dynamics of the cluster, as an important 
heating source is removed when the merged system leaves the GC. The detection of one IMRI would constitute a test 
of GR, as well as an irrefutable proof of the existence of IMBHs. 

Key words. (Galaxy;) globular clusters: general - gravitational waves - Methods: numerical - Stars: kinematics and 
dynamics 



1. Motivation 

IMBHs possibly exist in the centres of GCs vifith masses 
between M lO^"** Mq, if we assume that they follow 
the observed correlations be tween supermassive BHs and 
their stellar surroundings (see lMiller fc Colbertll2004l : lMilleij 
l2009l and references therein). Nonetheless, we do not have 
a clear detection of an IMBH, although there are favouring 
evidences that suggest that t hey should exist. In particu- 
lar, the work of iFarrell et"all ([2009, 20ll) derive a conser- 
vative lower limit of 500 Mq for an IMBH in the edge-on 
spiral galaxy ES0243-49 based on the variability of a X-ray 
source which has a maximum 0.2—10 keV luminosity of up 
to 1.1 X lO''^ ergs~^. Other recent interesting observational 
examples t hat point to the exi stence of these objects are 
the work of ISutton et al.l (j2012[ ). which evaluates a sample 
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of eight extreme luminosity ultra-luminous X-ray source 
candidates and state that the observed luminosities can be 
explained in t erms of IMBHs with masses in the range of 
10^ - 10'' Mr^. lNvland et al.l (|2012[ ) present Expanded Very 
Large Array observations of the dwarf lenticular galaxy 
NGC 404 and find that the most likely scenario is an accret- 
ing IMBH, ruling out other scenarios based on the observed 
X-ray and/or radio luminosities. 

Contrary to our Galactic Centre, we cannot use current 
observational techniques to detect a massive black hole in a 
GC. The radius of influence of an IMBH i s of a few arc sec- 
onds (|Peebleslll972l:IChaname et al.l[201(Tl:lMiller fc Colberti 
I2004D . For instance, for a 10'' Mq IMBH, the influence 
radius is of ~ 5" assuming a central velocity dispersion 
of a = 20kms~^ and a distance of ~ 5 kpc (see also 
iMiller &: Colbertll2004l for a similar example). 

There are other uncertainties, such as the availability of 
a bright reference star, which make the case rather improb- 
able. Hence, we would need the Very Large Telescope inter- 
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ferometer and one of the next-generation instruments, the 
VSI or GRAVITY (iGihessen et all 120061: lEisenhauer etaD 
[2008[) . In that case we could improve the astrometric accu- 
racy by an order of magnitude and thus wc would possibly 
be in the position of detecting the innermost kinematics of 
a GC around a potential IMBH. An interesting avenue to- 
wards the direct detection of an IMBH is GW astronomy. 
Additionally, these systems represent an excellent test of 
GR, since they are similar to extreme-mass ratio inspirals 
(jAmaro-Seoane et al.ll2007t) . In partic ular, space-borne de- 
tecto rs such as the ESA-led eLISA (jAmaro-Seoane et alj 
[20TI) or Chinese mission study options ( "ALIA" from now 



onwards, see iBender et al II2005I: ICrowder fc Cornishll2005t 



iGong et al.ll201ll) will be able to catch these systems with 
good signal-t o-noise ratios (SNR) if the GC is not further 
than z ~ 0.7 (jAmaro-Seoane et al.ll2Q12|) . 



2. Numerical tool and initial conditions 

Wc integrate the dynamical evolution of a globu lar cluster 
with Myriad (jKonstantinidis fc Kokkotas! |2Q10|) . a direct- 
summation A'^— body code that integrates all gravitational 
forces for all particles at every time step. The p rogramme 
uses the Hermite integration scheme (jAarse thl figgg. , 20031 
This requires computation of not only the accelerations, 
but also their time derivatives. Particles that are tightly 
bound or with very small separation are integrated usin g 
the time-symmetric Hermite scheme (iKokubo et al.l[T"998f) . 
which is a symplectic integrator that makes the numerical 
errors oscillate between two limits that can be controlled by 
the choice of the time step. The code uses post-Newtonian 
correcting terms to the Newtonian forces, including 1, 2 
and 2.5 order, as described for the first time in an A^— body 
code bv lKupi et al.l (I2006D (their equations 1, 2 and 3), as 
well as a recipe for gravitational recoil. The recoil velocity 
depends strongly on the mass ratio of the two holes, on 
the magnitude of their spins and on t heir directions with 
respect to the plane of the orbit (see e.g. lRezzo lial [2009l and 
references therein). The e quation that we have implemented 
in the code is taken from lLousto et ahl (l2010l ). 



V = {v,n + w± cos^) ei -I- v±_ sin -l- v\\ 63. (1) 

In the last equation, the indices _L and || stand for perpen- 
dicular and parallel directions with respect to the orbital 
angular momentum vector L of the binary, ei is a unit vec- 
tor and lies on the plane of the orbit connecting the two 
MBHs, with direction from the heavier to the lighter one. 
62 is also on the plane of the orbit, but perpendicular to 
ei, with direction such that ei, 62 and 63 construct an or- 
thonormal system, with 63 defined such that it is the unit 
vector parallel to i. ^ is the angle between the unequal 
contributions of mass and spin to the recoil velocity. We 
assign random, maximal spins to the stars of the GC, and 
in parti cular we initially gi ve the IMBH a spin a — S/A'P 
(see e.g. iLousto et al1l2010[ ) of 0.998. 

The number and masses of SBHs are of particular im- 
portance in our simulations. For this reason we created a 
large number of initial data for star clusters and evolved 
the s tars using the stellar evolution code sse (jHurlev et al.l 
l200Cl[). We used K roupa-like initial mass functions (IMF, 
see lKroupall200l with different values for the slopes ai 
(slope for m < 0.5Mq) and 012 (slope for m > O.SMq). We 



fixed the total number of stars to = 32768, the lower 
stellar mass limit to miow = 0.2Mq and the upper mass 
limit to TOuppor = 15OM0. Finally, we used values for the 
mctallicity Z ranging from 0.0001 to 0.02. For each model 
we created 100 sets of initial data and evolved the stars 
for 10 Myr, because this is the time in which all of the 
most massive stellar-mass BHs, the ones that could inter- 
act strongly and have an impact on the IMBH, form. In our 
models we assumed no supernova kicks, so that all SBHs 
formed are retained in the cluster. The number of SBHs 
depends strongly on the choice of the IMF and ranges from 
~ 20 to ^ 70 in our models. On the other hand the masses 
of the SBHs depend on the mctallicity and range from ^ 3 
Mq (for Z = 0.02) to - 27 Mq (for Z < 0.001). 

Since it is not possible to dynamically evolve all the ini- 
tial data sets we created, we picked 4 representative cases 
that produce low and high numbers of SBHs. We also picked 
a value Z = 0.001 for the metallicity as typical for a GC 
which resulted in the formation of SBHs with masses be- 
tween ~ 13 Mq and 27 Mq. In those models all stars with 
masses above 20 Mq have evolved off the main sequence at 
10 Myr. 

For our fiducial simulation A we choose slopes ai = 1.3 
and a2 ~ 2.4, which results in 62 stellar-mass BHs in the 
system. For the distribution of stars and BHs in the clus- 
ter we use a King profile (jKine Ji)66) with concentration 
parameter Wq ~ 7. The escape velocity at the centre of 
the clusters is ~ 17kms~^. No primordial mass segrega- 
tion is taken into account, so the BHs formed in all the dis- 
tances from the centre of the cluster. The absence of initial 
mass-segregation leads to an underestimate in the number 
of compact objects that will interact with the IMBH. The 
sizes of a star or a stellar remnant are taking into account 
only in the case of collisions. At the centre of the cluster we 
introduce an IMBH of mass M, = 500 Mq and correct the 
velocities of all stars and BHs of the GC to reach dynamical 
equilibrium. 

We run also three additional simulations with different 
parameters. Case B is like A but with M, = IjOOOMq 
and a2 = 2.5, which results in 52 SBHs; case C is like B 
but with a King parameter of 6 and 48 SBHs, and D is 
like A but with a King parameter of 6 and ai = 1.2 and 
Q!2 = 2.7, which result in 17 SBHs. Nevertheless, only case 
A had an IMRI; we will therefore in the remainder of the 
article focus on this case. In all simulations we turned off 
further stellar evolution, which, if turned-on, would form 
more SBHs in the GC. This is a simplification, which does 
dot have a significant impact on the dynamics of the IMBH 
and therefore on our results, since the low-mass SBHs that 
would form if stellar evolution was turned-on, would have a 
negligible infiuence on the IMBH, even if they merge with 
it. 



3. Dynamics of the system 

Initially, the IMBH interacts strongly with a sub-group of 
stars and BHs that contains approximately 20 members. As 
the system evolves, the members of this sub-group change. 
Soon, most of the stellar-mass BHs of the system sink to- 
wards the centre and start to interact with the IMBH and 
its environment. During this process, some of them receive 
big kicks due to 3-body interactions and are slingshot away 
from the centre of the cluster or GC itself. After T ^ 3 Myr 
the first stable IMBH-SBH binary forms. The companion of 
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the IMBH is a SBH with mass m,^ n = 23.9 Mq and the ini- 
tial semi-major axis of this binary is a ^ 88 AU. At T 9.2 
Myr this binary has a strong interaction with another SBH 
of the system. The interaction leads to a change of com- 
panion for the IMBH, which now builds a binary with a 
SBH of mass m,^ ig = 20.1 Mq. The initial semi-major axis 
of the new binary is a ^ 17.6 AU. This binary survives for 
nearly 40 Myr, but its characteristics vary significantly. At 
T ^ 49 Myr the semi-major axis changes to a ^ 5 AU, 
while the eccentricity increases to e = 0.965. At this point 
in the simulation, this binary interacts strongly with the 
second most massive SBH, which leads to a companion ex- 
change. The new binary has an initial semi-major axis of 
a ^ 6.55 AU and a very high eccentricity, of e = 0.999. 
The mass of the new companion SBH is to,. 2 = 26.54 M©. 
In Fig.([T|) we show the evolution of the semi- major axis 
and eccentricity for all of these binaries combined into a 
single curve. After some T ~ 13, 000 yr the binary merges 
and the resulting IMBH receives a random recoil velocity 
that depends on the mass ratio of the two members of the 
system and on the random spins that the code assigned 
to them. This "gravitational rocket" or recoil is such that 
the resulting velocity exceeds the escaping velocity and the 
merged system leaves the GC. This is due to the fact that 
we are using a low number of stars for the clusters; more 
realistic clusters will have larger escape velocities, so that 
the retained fraction of recoiling IMBH is larger and not 
well-represented by our case. We studied the distribution 
of recoil velocities for a merger of a binary similar to that 
of simulation A. We ran a two-body interaction lO'' times 
and calculated the recoil using equation [T] with different 
spin orientations and magnitudes for the two black holes. 
We found that the most probable recoil velocity for a bi- 
nary such as the one of case A peaks around 25 kms~^, 
with a probability of 21% that the merged system achieves 
velocities greater than 50kms~^, of the order of realistic 
escape velocities. 
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Fig. 1. Evolution of the semi- major axis and eccentric- 
ity of the different three binaries formed with the IMBH. 
Shortly after the beginning of the simulation, the IMBH 
builds a binary with the SBH with the 11th most massive 
mass, TO,^ 11. This corresponds to the first part of the curve 
(dashed blue curve). Later there is an interaction which 
leads to a companion exchange for the binary, the SBH 
with the 18th most massive mass, m.^ig. This binary lives 
for about 40 Myr. We can see that the two first binaries 
have phases of very high eccentricity, ebin ~ 1 ■, but not high 
enough to lead to a coalescence. The jumps in ebin indicate 
that the IMBH-SBH is still in a regime in which interac- 
tions with other stars play an important role. The system 
shrinks further and further until there is a three-body in- 
teraction. The binary is unbound and for a short period 
of time the IMBH has no companion, as indicated in the 
zoom-in subplots embedded in both, the upper and lower 
panels. Then the final binary forms, with the second most 
massive SBH. This binary is very hard and quickly losses 
energy via GWs radiation, which very efficiently leads to 
circularization and the final merger. 



In Fig. ^ we show the evolution of the distances of the 
10 most massive SBHs from the center. The SBHs inspiral 
the center very rapidly, as long as the IMBH exists in the 
cluster. Some of them escape the system, after passing very 
close to the central binary. After the IMBH merges with its 
binary companion SBH, the coalesced system leaves the GC 
and the trajectories of the remaining SBHs are not as steep, 
because they orbit the center of density of the GC without 
sinking rapidly into it. In Fig. ([3]) we show the Lagrange 
radii of the cluster during the simulation. We stop the sim- 
ulation at lOMyr after the ejection of the IMBH. The 
Lagrange radii increase instantly when the merged system 
leaves the system, but they start decreasing slowly shortly 
after that. This decrease will probably lead to a new core 
collapse of the system unless some hard binaries are formed. 
The IMBH-SBH binary represented a heating source in the 
center of the GC, transferring kinetic energy to the stars 
and SBHs that passed close from the density center. The 
lack of the binary makes the system core and internal radii 
decrease in size unless another heating source is formed and 
stops the shrinkage of the GC. Even though we integrated 
the system for another ^ lOMyr after the ejection, no such 
source formed. 



4. Gravitational waves from an IMRI 

In this section we follow the binary IMBH-SBH from the 
standpoint of emission of GWs. In Fig.Q we can see the 
evolution of the IMRI in a semi-major axis and orbital pe- 
riod - eccentricity plane. The binary enters the plot from 
the top with a high eccentricity, which places it very close 
to the innermost stable circular orbit, but the loss of energy 
quickly circularises it and drives it to lower eccentricities. 
As we discussed in the previous section, the binary forms 
with a very small initial semi-major axis, so that it hard- 
ens very efficiently. Hence, the binary follows closel y what 
we can expect from the approach of iPeter j ()1964|) . since 
the post-Newtonian terms lead the evolution of the sys- 
tem, which can be regarded as dynamically decoupled from 
the GC. It then enters the band of a LISA-like or ALIA 
detector with a significant eccentricity and the simulation 
is stopped when the semi- major axis is a = 5i?schw, the 
Schwarzschild radius of the IMBH. That is the moment at 
which the code assigns a recoil velocity to the merged sys- 
tem based on the spins of the two compact objects. 

In Fig.© we can see the same from the perspective of 
the characteristic amplitude and frequency of the waves. 
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Fig. 2. Distance R, to the density center of the GC of the 

ten most massive SBHs and the IMBH (sohd black hne). 
Strong interactions of the SBHs lead to ejection of four of 
them before the IMBH merges. They are removed from the 
simulation when R, > 10 pc and they are unbound with 
the GC. At r 47.7 Myr the IMRI leads to a coalescence 
that kicks the resulting merged system off the GC. 
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Fig. 3. Lagrangian radii showing the evolution of different 
mass fractions in the cluster: from the top to the bottom 
90, 80 ... 20, 10, 5, 3, 2, 1 and 0.1% of the total mass. The 
green rectangle shows the interval of time before and after 
the kick of the IMBH off the cluster. All mass curves suffer 
a jump at the moment of ejection. After the removal of the 
heating source from the center, the curves are flatter and 
their slopes start to decrease. 



We display the first harmonics in the appro ximation of 
Keplerian ellipses of iPeters fc Mathews! (jl963f ). 



5. Conclusions 

In this work we have investigated with a direct-summation 
code the evolution of GCs that harbours an IMBH in their 
center. The code uses relativistic corrections and a prescrip- 
tion for gravitational recoil. For one of the cases we find 
that an IMRI forms with a SBH due to close interactions, 
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Fig. 4. Inspiral of the SBH into the IMBH from the top to 
the bottom and from the left to the right in the eccentricity- 
semi-major axis plane. The red, solid curve starting at a 
very eccentric orbit shows the results of the Nbody simu- 
lation. The dashed, green region corresponds to the band 
of a LISA-like mission. Dashed, blue curves correspond to 
the trajectories due only to the emission of GWs in the 
iPeters fc Mathewi (|1963[ ) approximation. We also plot the 
corresponding merger timescales in the same approxima- 
tion in dashed, blue lines starting at 10^*^ years, and in 
solid, black lines the corresponding trajectories for evolu- 
tion by GW emission [Pct crs (196J) approximation. The 
black-shaded region on the right corresponds to the last 
stable circular orbit. Since the binary starts at a very high 
eccentricity, it basically follows one of the solid black lines, 
because it merges quickly and does not interact with other 
stars in the system. 



which leads to the ejection of the binary after coalescence. 
We follow the properties of the IMRI from a standpoint of 
GWs and the global evolution of the cluster. For z ^ 0.7 the 
energy loss in GWs is easily detectabl e by space-borne mis- 
sions such as a LISA-like observatory (jAmaro-Seoane et al.l 
120121 ) or ALIA in its 8 pc configuration. Moreover, the IMRI 
enters the bandwidth of the detectors with a very high ec- 
centricity, e = 0.9987, as with the EMRIs. One year be- 
fore the final coalescence, the system still retains a residual 
eccentricity of e '--^ 0.12. and ten minutes before merger 
of e ~ 10~^, which is detectable by data-analysis tech- 
niques (lAmaro-Seoane et all 120101 : iPorter fc Sesanal I2OIOI : 
iKev fc Cornishll2011h . 

IMRIs represent a test of GR, as well as a probe of 
space and time around massive black holes and also of the 
innermost kinematics of GCs to very large distances, of the 
order of a few Gpc. On the top of that, a successful detec- 
tion would represent very robust proof for the existence of 
IMBHs. The fact that the kick is making the merged system 
leave the GC is possibly an artifact of the low particle num- 
ber we used in the simulations, though in principle recoiling 
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Orb. evolution: Peters (64) 



GW amplitudes: Peters & Mathews (63) 
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Fig. 5. Characteristic amplitude he of the first harmon- 
ics of the quadrupolar gravitational radiation emitted dur- 
ing the inspiral of the IMRI. The numbers show the first 
four o f the h a rmon ics. The orbital evolution is calculated 
in the Peteri (11964 ) a pproxi mation and the amplitudes as 
in iPeters fc Mathew^ (Il963l ) . We assume the source is at 
a distance D = 1 Gpc. We indicate with a solid curve 
the noise curve \/ f Sh{f) for the ALIA detector with an 
armlength of 3 x 10^ m, a telescope diameter of 0.58 m, 
and a 1-way position noise of_8A/Hz pm; i.e. the 3H con- 
figuration of IXu. S.. et al.l (|201lf ). We also add the nois e 
curve for a LISA- like detector (in grey, iLarson et al.ll200d[) , 
with the Galactic bin ary white dwarf confusion background 
([Bender fc HilslllQQTl ). Note that the SNR is not given by 
the height above the curve, but by the area below it. For 
each panel we show the ratio i?p/i?s, the initial periapsis 
distance over the Schwarzschild radius of the system. We 
indicate the moments in the evolution for which the time 
to coalescence is 5, 1 yr, 1 hour, 10 minutes and 1 sec. 



velocities can be much higher than the escape v elocity of 
a clu s ter, of the order of ^ 50 kms ~^ (see e.g. iRezzo 
120091: iHoUev-Bockelmann et all I2008D . However, we have 
also demonstrated that there is a non-negligible statisti- 
cal probability that a similar case leads to a kick of the 
IMBH off a realistic GC. In spite of the code been ported 
to run on a PC with special-purpose hardware GRAPE, we 
can not cover a broader parameter space, nor study cases 
with a larger number of stars, or study the global dynami- 
cal evolution of the GC after the kick for a longer time. We 
plan on performing more accurate integrations on GPUs 
which will allow us to address the limitations we described 
above. This will allow us to investigate the potential global 
structure of the GC after the kick, since the impact on the 
cluster could in principle be a signature for the process. 
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